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We present results illustrating the effects of using explicit summation terms for the r~ 6 dispersion 
term on the interfacial properties of a Lennard- Jones fluid and SPC/E water. For the Lennard- Jones 
fluid, we find that the use of long-range summations, even with a short "crossover radius," yields 
results that are consistent with simulations using large cutoff radii. Simulations of SPC/E water 
demonstrate that the long-range dispersion forces are of secondary importance to the Coulombic 
forces. In both cases, we find that the ratio of box size Ln to crossover radius plays an important 
role in determining the magnitude of the long-range dispersion correction, although its effect is 
secondary when Coulombic interactions are also present. 



I. INTRODUCTION 

One of the fundamental problems in molecular sim- 
ulations is how to sum efficiently the pairwise inter- 
action potential u(r) that forms the basis for calcula- 
tion of nonbonded interactions. If the potential satisfies 
\u(r)\ < Ar~ m for m > 3, then the sum is guaranteed to 
be absolutely convergent; otherwise, it is at best condi- 
tionally convergent. Even in cases where the sum is abso- 
lutely convergent, however, explicit calculation of all the 
pairwise interactions, an 0(N 2 ) operation per time step, 
quickly becomes computationally unfeasible for molecu- 
lar dynamics simulations. Circumventing this difficulty 
usually involves one of two routes. In the first, a cutoff r c 
is introduced, outside of which interactions are neglected 
explicitly during the simulation. Long-range contribu- 
tions to chemical potential, surface tension, and other 
thermodynamic properties are calculated by integrating 
from the cutoff radius r c to infinity, assuming the radial 
distribution function g(r) equal to unity. Certain cases, 
however, require large cutoff radii, which leads to compu- 
tational inefficiency because of the large number of pairs 
needed to calculate the total pairwise sum. 

An alternative method involves transformation of the 
pairwise interactions beyond some "crossover" distance 
(which we designate as r^) into a more readily-solved 
problem. Although truncation of the pairwise potential is 
more computationally efficient, and has in fact been used 
in parametrization of the SPC/E,ETIP3PPand TIP4P^ 
models of water, excluding interactions from considera- 
tion in the simulation can introduce significant errors in 
the values of thermodynamic properties obtained from a 
simulation. In such cases one generally relics on Ewald 
summations, a technique which replaces the long-range 
calculations by a sum of replicas of the central box. We 
consider the effects of including explicit summation of 
long-range dispersion forces. These long-range forces are 
normally neglected in one-phase systems, but become of 
importance in multiphase systems or when studying in- 
terface phenomena. 

Following earlier work by Williams 3 and Perram and 
co-workersP Karasawa and Goddard used long-range 
summations of dispersion forces to derive thermodynamic 
properties of argon and of NaCl crystalsPMore recently, 



several other investigations into the use of long-range 
summation of dispersion forces have been reported. Ou- 
Yang et al. have used long-range summation to study 
the phase equilibrium of Lennard- Jones fluids P Lopez- 
Lemus and co-workers have used Ewald summation to 
study the density of n-alkanes!^ Shi et al. have extended 
the particle-particle particle-mesh (PPPM) implementa- 
tion of Hockney and EastwoocP to include Ewald sums 
for the dispersion term, and used it to study the phase en- 
velope of water PEssmann and co-workers have created a 
similar mesh-based formulation for Ewald summations 
Greater emphasis has been placed on how to use long- 
range corrections, as applied to, for example, constant- 
pressure^ and inhomogeneous^ simulations. Deserno 
and Holm considered the effects of various parameters 
involved in performing the Ewald sum matio n on its accu- 
racy for several mesh implementations}*^ 3 -^ while Pollock 
and Glosli have evaluated the efficiency of both tradi- 
tional Ewald summations and mesh implementations.^ 



In this paper, we present results on the use of Ewald 
summations to handle the r -6 dispersion term. We first 
illustrate the effect of incorporating the long-range sum 
by studying a model Lcnnard-Jones fluid; we report the 
liquid- vapor coexistence data, as well as the surface ten- 
sion. Next, we demonstrate the effects of long-range sum- 
mation of the dispersion term in ionic systems by simulat- 
ing SPC/E water. In both cases, we compare the Ewald 
simulation results to results of simulations with trun- 
cated Lennard- Jones potentials. We also examine the 
role long-range corrections of the surface tension plays in 
both types of simulations. 



In Section [IT] we briefly outline the equations that must 
be implemented to evaluate the Ewald sum for potentials 
of the form /(r) ~ Ar~ 6 . We then demonstrate its appli- 
cation to a Lennard-Jones fluid and to water, which has 
both Coulombic and dispersion Ewald sums, in Sections 



III and IV We consider issues related to the performance 
of the various algorithms for performing long-range sum- 
mations in Section [V] Finally, we offer our conclusions in 
Section rVTl 
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II. MATHEMATICAL FORMULATION OF THE 
EWALD DISPERSION SUM 

Incorporation of long-range summation of LJ poten- 
tials, 



Uij(r) = 4e; 



D 



(1) 



requires treatment of the second, attractive term; the 
repulsive r -12 term can be neglected, particularly as r 
becomes large. We now briefly summarize the equations 
governing long-range summation of the attractive r~ 6 or 
dispersion portion of the Lennard- Jones potential. Note 
that the following also applies to any other potential with 
attractive r~ e contributions, such as a Buckingham exp- 
6 potential, provided that excluded volume interactions 
fall off fast enough. 

Karasawa and GoddarcP' derive expressions for the dis- 
persion part of energy, force, and stress using Ewald sum- 
mations. For the energy they find 
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where V represents the volume of the simulation cell, h 
the reciprocal lattice vector, h — |h|, a = |r.; — Yj~ Rl|/?7, 
and b — hri/2. The first contribution to the energy sums 
over lattices L and particles i and j. The parameter R^ 
represents the lattice translation vector. The parameter 
r\ represents the range of interactions handled in recip- 
rocal space; for a given value of r, the larger the value 
of 77, the more interactions that are included in recipro- 
cal space. We define the quantity as the "crossover" 



which represents the maximum distance for which inter- 
actions are summed using the first term in Eq. [2] Note 
that in combination with the number of reciprocal lat- 
tice vectors n^ ec determine the numerical accuracy with 
which the lattice-sum Hamiltonian is evaluated. Using a 
finite n^ ec in combination with a finite results in an 
approximate solution for the lattice sum. Conversely, the 
solution is exact when either one is infinite. The excluded 
volume uses the same cutoff as the crossover radius, so 
no additional cutoff for the r~ 12 term is needed. The 
structure factor S^h), a complex number, is defined as 



(3) 



assuming geometric mixing of dispersion constants 
(Bij = ^BaBjj = bibj). This type of mixing applies 
to OPLS potentials by Jorgensen et aiP^ On the other 
hand, potentials like CHARMlvf^ use arithmetic mixing 
for sigma cross-terms (ov,- = \(Tu+ ^Pjj), in which case 
it is necessary to expand the dispersion constant, 
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with ft^fc = an d Cfe = (,). This redefines the 

structure factor product as 



S 6 (h)S 6 (-h) = 2^S 6 , fc (h)S 6 , 6 _ fc (-h), (5) 
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with 



Se,fe(h) = ^6j,fcexp(-zh • r,). 
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Note that symmetry requires the calculation of only four 
out of the seven structure factors products in Equation 
[5] The force on atom k is given by 



f 6 ,fe = -^^BkjiTk-rj -R L )(6a- 8 + 6a- 6 + 3a- 4 + a- 2 ) 
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The instantaneous stress is given by 
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III. METHODOLOGY 



tween 320 and 1404 k-vectors. 



A. Lennard-Jones fluid 

All molecular dynamics simulations were performed 
using the LAMMPS simulation package^ with added 
Ewald dispersion summations. All Lennard-Jones pa- 
rameters can be expressed in terms of the distance tr^ = 
a, the root of the Lennard-Jones potential Q; ey, the 
well-depth of the potential; and the time, r = a(raj 'e) 1 / 2 , 
where m is the mass of the atom. Initial configurations 
were created by randomly placing either N part i C ies — 4000 
or 16000 particles in simulation boxes with respective 
geometries of L x = L y = L\\ = ll.Olcr or 22.02er, and 
L z = L± = 44.04cr. Avoidance of system size effects due 
to large cutoffs motivated the choice of using system sizes 
of either 4000 or 16000 particles. 

Thus created, configurations were equilibrated for 500r 
at T* = k B T/e = 0.7 with At = 0.005r, using three- 
dimensional periodic boundary conditions. Once equili- 
brated, liquid-vapor interfaces were created by changing 
the size of the box in the z-direction to L± = 176.16<r 
without displacing the particles. Once strained, the sys- 
tem was equilibrated for 500r at T* — 0.7. Subsequent 
higher temperatures used this state as initial configura- 
tion. To reach the desired temperature, the system was 
slowly heated over a period of 500r; the system was then 
allowed to equilibrate at the final temperature for 250r. 
Figure [T] shows a typical snapshot of a 16000 atom sys- 
tem. 

Once at the desired temperature, total system pres- 
sure tensors and z-directional density profiles were mea- 
sured every 5r for 1250 to 5000r. All simulations used 
a Nose-Hoover thermostat with relaxation time 10 r. 
Studied temperatures were T* G {0.70,0.85,1.10,1.20}. 
Values for r\ follow from the choice of precision through 
rj = (1.35 — 0.15 In P)/r*, where P is the precision. 1 - 
Ewald dispersion summations used a precision of 0.05. 
Changing precisions to 0.01 did not change the results. 
For a given system size, the number of k-vectors em- 
ployed in each simulation decreases with increasing r^; 
the Lennard-Jones simulations described here use be- 



B. SPC/E water 

The dominant contribution to the surface tension 
of water is from the electrostatic interactions, which 
are much stronger than the Lennard-Jones interactions. 
However, the traditional Ewald method or the PPPM 
technique of Hockney and EastwoocP is generally used 
for the long-range electrostatics only. Long-range cor- 
rections are applied to the surface tension 19 ! 2 ™ a nd other 
system properties^ to account for the effects of truncat- 
ing the Lennard-Jones potential. 

The SPC/E water modeP was used for all simulations 
described below, as we found in a previous studj^ that 
the SPC/E model showed acceptable agreement with ex- 
perimental data while avoiding the computational ex- 
pense associated with four-point water models. The 
SHAKE technique^ was used to constrain the bond 
lengths and bond angles. The equations of motion were 
integrated using the Verlet algorithm; a Nose-Hoover 
thermostat with a relaxation time of 500 fs was used to 
control the temperature during the simulations. All sim- 
ulations used a timestep of At = 1.0 fs. 

To determine the effect of the long-range Ewald sum- 
mation of the dispersion term on the liquid- vapor equilib- 
rium properties of water, water molecules were randomly 
placed in a periodic, rectangular box. In all cases, the di- 
mensions of the box in the directions parallel to the inter- 
face were equal: L x = L y = Ln, where Lu varied between 
3.65 nm and 7.74 nm. For all simulations, the dimension 
perpendicular to the interface was L z = Lj_ = 5.0 nm. 
The initial density in each box was p — 1.00 g cm~ 3 . 

For each different set of Lennard-Jones and electro- 
static cutoffs, systems were allowed to equilibrate at 300 
K or 400 K for 500 ps. The resulting configuration was 
then placed in a box with L z = L± = 15.0 nm, and 
allowed to equilibrate for 1 ns to establish the liquid- 
vapor interface. Production runs of at least 1 ns then 
followed, with pressures measured every timestep and po- 
sitions stored every 5 ps. 
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FIG. 1: Snapshot for a 16000 atom monomeric Lennard- Jones fluid at coexistence with cutoff r c = 7.5<r at T* — 1.10 shown in 
the xz-plane. 



With one exception, for the systems which included 
long-range summation of the Lennard- Jones dispersion 
term, the Ewald summation was also used for the electro- 
static interactions. In such cases, the crossover radius for 
both the LJ and electrostatic interactions were identical. 
For the systems in which the Lennard-Jones potential 
was truncated, electrostatic interactions were calculated 
using the PPPM methodPThe PPPM method yields re- 
sults that are indistinguishable from Ewald summations 
(within simulation uncertainty), but have much faster 
running times for large systems. The root-mean-squared 
precision of the algorithm was on the order of 10~ 3 . We 
showed in our previous paper that the interfacial proper- 
ties of SPC/E water do not show a significant dependence 
on the accuracy of the k-space mesh, provided that the 
system is allowed to equilibrate sufficiently In our previ- 
ous studies of the surface tension of water models,^ we 
found that k-space meshes of 5 x 5 x 20 through 20 x 20 x 80 
yielded equivalent results. The PPPM simulations re- 
ported here used meshes of size 10 x 10 x 27; the Ewald 
summations used a total of 2046 vectors. 



The outer factor of 1 /2 in Eq. ( 10 1 accounts for the pres- 
ence of two liquid- vapor interfaces. 

As discussed above, the surface tension measured 
by simulations, 7 S i m , tends to be lower than the 
experimentally-measured surface tension. Improved 
agreement with experimental data by incorporating the 
tail correction, 7t a ii: 

7 = Ip+ltaiU 

where 7t OJ ; can be determined fronii2122l 
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In Eq. 11 U (r) is the pairwise potential, g (r) is the ra- 



dial distribution function, p (z) is the observed interfacial 
profile, and pa ( z ) is a Gibbs dividing surface: 
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C. Surface tension and tail corrections 

There are two primary methods used to compute the 
surface tension using molecular simulations. One ap- 
proach, developed by TolmaiPS and refined by Kirkwood 
and Buff]2£l computes the surface tension as an integral 
of the difference between the normal and tangential pres- 
sures p± (z) and pii (z): 



Eq. [TT] represents the contribution to the stress from the 
region r > r c which is not captured by simulation; the 
term containing the product of densities is an approx- 
imation for the true pair density p^ required by the 
Kirkwood-Buff theory.^ Following our earlier work on 
tail corrections in the surface tension of water models, 
we fit the density profile to an error function profile to 
simplify the calculations.^ 
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where, in our geometry, p±(z) — p z (z), and p\\{z) — 
(p x (z) + p y (z))/2. Away from an interface, p± = pn and 
the integrand vanishes. Therefore the nonzero contribu- 
tions to the integral in Eq. ^ come from the interfa- 
cial region. In the case of an interface between two fluid 
phases, the integral in Eq. Q can be replaced with an 
ensemble average of the difference between the normal 
and tangential pressures: 



iv = y (p±-P\\) = y 
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IV. RESULTS AND DISCUSSION 
A. Lennard-Jones fluid 

Surface tensions at interfaces in two-phase systems 
are sensitive to long-range interactions, even when sim- 
ple Lennard-Jones interactions are concerned. This phe- 
nomenon becomes especially pronounced close to the crit- 
ical point, where large fluctuations in density dominate. 

The interfacial properties of the various L J simulations 
are collected in Table [T] The effects of using a short cut- 
off are immediately apparent from this table and Figure 
[2] at all temperatures studied, there is a significant dif- 
ference between the calculated densities p* , pressure p* , 
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FIG. 2: Temperature dependence of the surface tension for 
a monomeric Lennard- Jones fluid at coexistence with cutoff 
r c — 2.5a (square), 5.0a (triangle), 7.5a (diamond) andEwald 
summation crossover = 5. Oct (star; overlaps with cutoff 
r c = 7.5a data). Lines are added as a guide to the eye. Re- 
ported surface tensions calculated by means of using a cutoff 
include tail corrections. 



and surface tension 7* = 7a 2 /e observed for a cutoff of 
r c = 2.5ct and cutoffs of r c = 5. Oct or r c = 7.5a. The 
liquid density of the r c = 2.5a simulation is statistically 
lower than the liquid densities at larger cutoffs; the va- 
por density is likewise statistically higher. The relative 
disagreement increases in size with increasing tempera- 
ture, as can be seen in comparing the density plots of the 
r c = 2.5a and r c = 7.5a simulations for T* = 0.70 and 
T* = 1.10, as shown in Figure [3] The interfacial region 
is much broader for the smaller cutoff, indicating that 
the system is near its critical temperature T* lit . Litera- 
ture confirms this observation: the critical temperature 
of Lennard- Jones fluids is approximately T c * rit ss 1.18S 26 
for r c — 2.5a, which represents about an 11 percent 
decrease in the critical temperature when compared to 
T c * rit « 1.316 for r c approaching innnityj23 

Measured densities for Ewald dispersion simulations 
vary by not more than 2 percent for the liquid phase and 9 
percent for the gas phase. Furthermore, good agreement 
is observed between densities of simulations using cutoffs 
of r c = 7.5a and corresponding densities obtained from 
Ewald dispersion simulations with — 5. Oct. Similar 
trends exist in the surface tension 7*. Simulations using 
a cutoff yield surface tensions that at each temperature 
increase with the cutoff r c . As the temperature increases, 
the observed values of 7* tend towards zero, which is 
another indicator that truncating the cutoff induces a 
lower critical temperature. 

For simulations using Ewald summations, differences 
due to changes in crossover radii are much smaller. In 
fact, for T* — 0.85 and T* = 1.10, there is essentially no 
difference between the values of 7* for crossover radii of 
4.0ct and 5. Oct; the results are identical within the uncer- 
tainty of the calculation. 

Mecke et alP^ confirm our results when implicitly 
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FIG. 3: Density profiles of a monomeric Lennard- Jones fluid 
at coexistence with cutoff r c — 2.5a (dashed) and Ewald sum- 
mation crossover = 4. Oct (solid) at a) T* = 0.70 and b) 
T* = 1.10. 



adding self-adjusting localized long-range corrections 
during simulation using on-the-fiy density distributions 
instead of long-range correction through Ewald summa- 
tions. Their method requires a priori knowledge of the 
spatial distribution of the phases and breaks down for 
three dimensional phase distributions. Application of 
Ewald dispersion summations circumvents this require- 
ment and is applicable to any type of phase distribution. 



B. SPC/E Water 

Our results for the water simulations are collected in 
Table [Til The most obvious difference between the results 
for Lennard- Jones fluids and those obtained for water is 
the difference in the relatively small importance of adding 
the long-range Lennard- Jones corrections to the potential 
for water. 

We examined a number of different systems to analyze 
the effects of the r~ 6 dispersion term, comparing systems 
with a cutoff to systems using Ewald summation for the 
dispersion sum. We first examined systems with uniform 
size and electrostatic crossover = 10 A; no clear trend 
was evident. Keeping L x fixed while increasing both cut- 
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TABLE I: Interfacial properties of Lennard- Jones fluid a 



Method ^particles 


v c j a 




T* 




J°vap 




7<* -i 
'tail 


7* 


Cutoff 16000 


2.5 




0.700(0) 


0.7868(1) 


0.0079(0) 


0.557(17) 


0.252 


0.809 




5.0 




0.700(0) 


0.8355(1) 


0.0026(0) 


1.033(17) 


0.104 


1.137 




7.5 




0.700(3) 


0.8401(1) 


0.0019(0) 


1.088(14) 


0.052 


1.140 




2.5 




0.850(0) 


0.7019(1) 


0.0274(0) 


0.331(05) 


0.183 


0.516 




5.0 




0.850(0) 


0.7675(1) 


0.0103(0) 


0.713(14) 


0.084 


0.797 




7.5 





0.850(1) 


0.7732(1) 


0.0100(0) 


0.762(08) 


0.043 


0.805 




2.5 




1.100(0) 


0.4780(3) 


0.0847(2) 


0.026(11) 


0.016 


0.042 




5.0 




1.100(0) 


0.6383(1) 


0.0453(1) 


0.259(12) 


0.050 


0.309 




7.5 




1.100(0) 


0.6391(1) 


0.0514(1) 


0.304(08) 


0.028 


0.332 
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3.0 


0.850(5) 


0.7663(1) 


0.0076(1) 


0.778(44) 




0.778 






4.0 


0.851(5) 


0.7722(1) 


0.0078(1) 


0.817(32) 




0.817 






5.0 


0.850(4) 


0.7742(1) 


0.0072(1) 


0.828(49) 




0.828 






3.0 


1.100(1) 


0.6353(2) 


0.0460(1) 


0.333(20) 




0.333 






4.0 


1.100(1) 


0.6346(2) 


0.0413(1) 


0.373(38) 




0.373 






5.0 


1.100(1) 


0.6452(2) 


0.0429(1) 


0.366(46) 




0.366 


16000 




5.0 


0.700(0) 


0.8414(1) 


0.0020(0) 


1.141(11) 




1.141 






5.0 


0.850(0) 


0.7750(1) 


0.0098(0) 


0.818(15) 




0.818 






5.0 


1.100(0) 


0.6383(1) 


0.0553(0) 


0.302(11) 




0.302 






5.0 


1.200(0) 


0.5627(2) 


0.0942(0) 


0.156(14) 




0.156 



^Digits between parentheses indicate the uncertainty in the last significant figure. 
b Above the critical temperature. 



offs uniformly likewise yielded no systematic behavior. 
Increasing r c and rjf while keeping the ratio L x /r c con- 
stant shows the expected increase in surface tension 7 s ; m 
with increasing cutoff. For the corresponding Ewald sim- 
ulations, the surface tension 7 S j m shows a much smaller 
variance with respect to changes in crossover and sys- 
tem size. Examining the overall surface tension, which 
incorporates the tail correction 7t a ;i, we find that there 
is significantly greater variations in the simulations in- 
corporating tail corrections for the r~ e term than those 
using Ewald summations for the long-range interactions. 

At 300 K, Ewald summation of the dispersion interac- 
tion yields a significantly improved estimate of the liquid- 
phase density pu q , particularly as the crossover distance 

increases. By contrast, at 400 K, even for a cutoff of 
r c = = 16 A, the disagreement between the simula- 
tion and experimental densities is approximately 3 per- 
cent, which reflects the disagreement of the commonly 
used nonpolarizable water models with e xperim ental re- 
sults at temperatures well above 300 KPSBSl However, 
differences between experimental and simulation results 
are significantly smaller when using Ewald summation 
than when employing a cutoff. 

Comparing the total surface tensions 7 to the experi- 
mental results, however, we find that there is still a sig- 
nificant discrepancy between the two sets of measure- 
ments. Although the inclusion of the long-range correc- 
tions yields an improved estimate compared to our pre- 
vious work comparing the surface tensions of three- and 
four-point water models,^ the simulated surface tensions 
are still underestimated by approximately 10 to 15 per- 
cent. These results suggest that the SPC/E water model 
is incapable of accurately modeling interfacial properties, 



particularly at temperatures away from 300 K (the tem- 
perature at which the model was parametrized). 

The relatively poor performance of the SPC/E water 
model may be an artifact of its original derivation, as 
the model is designed to perform well as bulk liquid, but 
is known to perform poorly when considering the vapor 
phase as a result of its dipolc moment, which is signifi- 
cantly larger than its experimental value. Such problems 
may be inherent to most non-polarizable models, which 
include most of the commonly used three- and four-point 
modelsp] polarizable models should perform better due 
to their adaptive nature. However, our choice of the 
SPC/E water model was driven by its simplicity and the 
desire to illustrate the effects of a charged system on the 
performance of long-range dispersion summations. 



V. ALGORITHMIC PERFORMANCE 

While Ewald summations are the standard method for 
handling long-range forces in Monte Carlo simulations, 
for molecular dynamics simulations, Ewald summations 
are practical only for a relatively small numbers of par- 
ticles. Pollock and Gloslp^have shown that for systems 
containing roughly 1000 particles or more, PPPM sim- 
ulations outperform simulations that use explicit Ewald 
summation. To reduce computational complexity, molec- 
ular dynamics simulations typically implement Fourier 
transform methods using meshes, such as the particle- 
particle particle-mesh (PPPM) method introduced by 
Hockney and EastwoodP The PPPM method replaces 
the exact evaluation of the long-range interactions with 



7 



TABLE II: Interfacial properties of SPC/E water 



T 


Method 






N w 






(Pvap) 


7sim 


Ttail 


7 


K a 




(A) 


(A) 




(A) 


(g/cm 3 ) b 


(g/cm 3 ) 


(mN/m) b 


(inN/m)' 


(mN/m) 


300 


Exptl. d 


— 


— 






0.9965 


3 x 10~ 5 


— 


— 


71.7 




PPPM 


10 


10 


2200 


36.5 


0.9889 


2 x 10" 5 


55.2 


6.4 


61.6 






12 


10 


2200 


36.5 


0.9924 


0.0000 


59.1 


4.6 


63.7 






12 


12 


2200 


36.5 


0.9909 


3 x 10 -5 


59.7 


4.6 


64.3 






12 


12 


3200 


43.8 


0.9917 


4 x 10~ 5 


57.1 


4.6 


61.7 






16 


10 


2200 


36.5 


0.9954 


0.0000 


58.4 


2.7 


61.1 






16 


16 


2200 


36.5 


0.9950 


3 x 10 -5 


56.1 


2.7 


58.8 






16 


16 


5600 


57.9 


0.9940 


0.0001 


59.4 


2.7 


62.1 






20 


20 


10000 


77.4 


0.9954 


3 x 10 -5 


59.6 


1.7 


61.3 




Ewald 




10 


2200 


36.5 


0.9954 


0.0006 


61.8 




61.8 








12 


2200 


36.5 


QQ56 


4 x 10 — 5 


61 6 




61.6 








12 


3200 


36.5 


9957 


o ooon 


61 1 




61.1 






12'' 


12 


3200 


43.8 


0.9912 


0003 


56.5 


4.6 


61.1 








16 


2200 


36.5 


QQ66 


0000 


59 6 




59.6 








16 


5600 


57.9 


QQ58 


0001 


60 5 




60 5 


400 


a 

Exptl. 










0.9374 


0.0014 






53.6 




PPPM 


10 


10 


2200 


36.5 


0.9076 


0.0012 


41.1 


5.0 


46.1 






12 


12 


3200 


43.8 


0.9106 


0.0014 


42.1 


3.7 


45.8 






12 


12 


2200 


36.5 


0.9106 


0.0015 


42.5 


3.7 


46.2 






16 


16 


2200 


36.5 


0.9138 


0.0012 


43.5 


2.2 


45.7 






16 


16 


5600 


57.9 


0.9143 


0.0012 


40.9 


2.2 


43.1 




Ewald 




10 


2200 


36.5 


0.9154 


0.0014 


43.8 




43.8 








12 


2200 


36.5 


0.9154 


0.0014 


43.6 




43.6 








16 


2200 


36.5 


0.9160 


0.0009 


43.0 




43.0 








16 


5600 


57.9 


0.9165 


0.0010 


44.1 




44.1 



a Average temperatures agree with the thermostat temperature within 0.02 K. 

Uncertainties for densities are less than 0.002 g/cm 3 ; for surface tensions, between 2.0 and 3.0 mN m 
c Estimatcd using tanh density profile. 
d Data taken from Ref. 1311 

c Long-range summation performed only for Coulombic interaction. 



an approximate treatment through the use of the fast 
Fourier transform over a discretized mesh. This change 
reduces the complexity of computing the forces from 
0(N 2 ) to 0(N log N) per time stepPl 

Besides the number of particles, the computational 
performance of Ewald summations also depends on the 
number of k-space vectors, which is derived from the cho- 
sen precision and crossover radius. We compare perfor- 
mance of Ewald summations to cutoff simulations in two 
ways, using timing data reported in Table Compar- 
isons between two simulations with equal number of par- 
ticles and equal cutoff and crossover radius shows that 
Ewald summations are typically 2^ to 3 times slower. 
However, the primary objective is to obtain results which 
approximate those that can be obtained without us- 
ing a cutoff: we find simulations of 4000 particles with 
a crossover radius of 5. Oct have the same accuracy as 
simulations of 16000 particles with a cutoff radius of 
r c = 7.5ct. The timing comparison in this case shows 
that Ewald summations are about l| times faster. Ewald 
summations have the added benefit of negating the need 
for post-processing of tail corrections. This leads us 
to conclude that, regarding performance and accuracy, 
Ewald dispersion summations of Lennard- Jones simula- 
tions are preferable to using cutoffs combined with tail 
corrections. 

Although for smaller systems, Ewald summations are 
more efficient than mesh methods, as system size in- 
creases, which method is preferred can depend upon 



architecture- and problem-specific factors. This becomes 
even more evident when one considers problems also in- 
cluding Coulombic interactions, as is the case with water, 
as shown in Table |IV| which shows running times for a 
system of 6600 atoms on four processors. The fastest 
running time for the Ewald simulations (357.0 ms) is sig- 
nificantly slower than the slowest running time for the 
PPPM simulation (308.6 ms), which uses a substantially 
larger cutoff radius (16 A versus 12 A). It should be 
noted, however, that only temperatures well below the 
critical temperature were studied. In general, we find 
that the PPPM simulation running time is dependent 
on the larger of the LJ cutoff and the crossover radius, 
while selection of an optimal Ewald crossover is compli- 
cated by the existence of an optimal crossover radius. 
The slower execution of the Ewald algorithm, however, 
does not lead to a major improvement in accuracy: the 
difference between the surface tensions for the PPPM al- 
gorithm with tail correction and the Ewald algorithm for 
the various systems is roughly 1 percent. Incorporation 
of the long-range dispersion interactions into a PPPM or 
other mesh-based formulation will be desirable for large 
systems. 



VI. CONCLUSIONS 

Simulated data shows the importance of long-range ef- 
fects on pressure, density, and surface tension at all tern- 
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TABLE III: Running times for LJ simulations with cutoffs versus Ewald summations 



Time per timcstcp (ms) 

Npar-ticles = 4000 N par . tic u s = 16000 



T* 


Method 


rc (o-) 


(4 procs) 


(8 procs) 


0.70 


Cutoff 


2.5 


5.2 


18.4 




Cutoff 


5.0 


30.4 


40.6 




Cutoff 


7.5 


98.7 


135.0 




Ewald 


5.0 


95.1 


124.6 


0.85 


Cutoff 


2.5 


5.4 


6.7 




Cutoff 


5.0 


33.1 


37.8 




Cutoff 


7.5 


115.1 


129.0 




Ewald 


5.0 


87.7 


121.9 


1.10 


Cutoff 


2.5 


3.5 


6.4 




Cutoff 


5.0 


25.0 


35.9 




Cutoff 


7.5 


76.2 


119.7 




Ewald 


5.0 


67.3 


103.5 



TABLE IV: Running times (ms per timestep) for water simulations at 300 K using PPPM and Ewald summation for 6600 
atoms on 4 processors 







Crossover radius 




Method 




r* = 10 A 


r* = 12 A 




= 16 A 


PPPM, r c = 


10 A 


93.3 


143.8 




289.4 


PPPM, r c = 


12 A 


140.9 


144.8 




296.4 


PPPM, r c = 


16 A 


277.7 


292.4 




308.6 


Ewald, r c = 




459.0 


357.0 




423.5 



peratures, but especially close to the critical point. For 
simulations of a homogeneous bulk system, cutting off 
the dispersion term at r c is sufficient provided that the 
pair correlation function g(r) for r > r c is close to 1.0. 
In this case, the forces for distances r > r c cancel, and 
the contributions to the pressure and energy for r > r c 
can be easily accounted. However, for two-phase systems 
such as those studided here, cutting off the dispersion 
contribution to the pair potential have strong effects on 
the properties of the system. Part of the missing contri- 
bution to the surface tension can be included by using 
tail corrections, but since truncating the potential at r c 
can have a strong effect on the liquid and vapor densities, 
particularly near the critical point, inclusion of the tail 
correction for the surface tension only accounts for part 
of the missing contribution. This supports the notion 
that description of long-range effects by means of Ewald 



summation is preferential close to the critical point. The 
long-range effects become less pronounced when further 
removed from the critical point. Comparison between the 
Lennard- Jones fluid and water results show that the im- 
portance of long-range dispersion effects become less pro- 
nounced when stronger contributions, such as Coulombic 
terms, are included. In this case, the latter is the dom- 
inant contribution, while the former plays a secondary 
role. 
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